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ABSTRACT 

We investigate the formation process of planetesimals from the dust layer by the 
gravitational instability in the gas disk using local iV-body simulations. The gas is 
modeled as a background laminar flow. We study the formation process of planetesimals 
and its dependence on the strength of the gas drag. Our simulation results show that the 
formation process is divided into three stages qualitatively: the formation of wake-like 
density structures, the creation of planetesimal seeds, and their collisional growth. The 
linear analysis of the dissipative gravitational instability shows that the dust layer is 
secularly unstable although Toomre's Q value is larger than unity. However, in the initial 
stage, the growth time of the gravitational instability is longer than that of the dust 
sedimentation and the decrease in the velocity dispersion. Thus, the velocity dispersion 
decreases and the disk shrinks vertically. As the velocity dispersion becomes sufficiently 
small, the gravitational instability finally becomes dominant. Then wake- like density 
structures are formed by the gravitational instability. These structures fragment into 
planetesimal seeds. The seeds grow rapidly owing to mutual collisions. 

Subject headings: gravitation, instabilities, methods :n-body simulations, planets and 
satellites:formation 

1. Introduction 

In the standard scenario of planet formation, planetesimals are the precursors of planets. 
Their formation process is one of the unsolved problems of the planet formation theory. Beginning 
with micron-siz ed dust grains, they grow to centimeter size in a protoplanetary disk via collisional 



agglomeration (|WeidenschillingHl980l : iBlum &: Wurmll2000l ). The least understood growth phase is 
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that of growth from centimeter size to kilometer size. Since the gas drag of meter-sized aggregates 
is weak, they are slightly decoupled from gas with the sub-Keplerian velocity. The resulting head- 
wind decreases their angular momentum and causes inw ard drift, whose timescale is about a few 
hundred years (jAdachi et al.lll976l ; IWeidenschillina 119771 ) . The growth process in this stage must 
be faster than this drift timescale. Another problem is the sticking probability of the aggregates. 
In this stage, sin ce the threshold velocity for the d estruction is low, their collisions may lead to 



destruction (e.g.. iBlum &; Wurm 



2000 



Sironol 12004 ). 



A very thin and dense layer of settled dust aggregates in the mid-plane of the protoplane- 
tary disk is gravitationally unstable. Then the gra vitational collapse of the dust layer occu rs, and 
kilometer-sized planetesimals are formed directly (|Safronovl Il969i ; iGoldreich &: Ward! Il973l ) . This 
scenario has the advantage of a very rapid formation timescale, which is on the order of the Ke- 
plerian time, thus avoiding the migration of meter-sized aggregates. However, turbulence in the 
protoplanetary disk may prevent dust from settling to the mid-plane. As dust settles in the proto- 
planetary disk, a vertical shear velocity develops because the dust-rich gas in the mid-plane rotates 
with the velocity closer to the Keplerian velocity than that of the dust-depleted gas, which may 
cause turbulence. The turbulence can mix dust with gas and prevent the d ust layer from settling 
into a dense enough sheet for gravitational instability (| Weidenschillinei 1 1980 ) . Problems concerning 
the gravitational instability of the dust layer still remain unsolved. 

There are two methods for calculating the dynamics of dust particles in gas: fluid simulation 
and A^-body simulation. One can consider the particles a s virtually fluid and trea t them with the 
Euler equation or Navier-Stokes equation. For example, lYamoto Sekiyal (120061 ) performed the 
two-dimensional numerical simulation to investigate the density evolution of the dust layer due to 
gravitational instability. They assumed that the dust layer is axisymmetric with respect to the 
rotational axis. They treated the dust component as a pressure-free fluid because the velocity 
dispersion of dust is negligible if the stopping time i s t op is sufficiently short, where the stopping 
time t st0 p is the characteristic timescale for a particle to stop in gas. They foun d that the dust 
layer becomes extremely thin if the stopping time is long. IWakita &: Sekiyal (120081 ) performed the 
numerical simulation of the gravitational instability of a two-dimensional thin disk and investigated 
the non- axisymmetric modes. 

On the other hand, we can investigate the dynamics of particles in gas by using A^-body simula- 



tions (jTanga et all2004l ; lMichikoshi et al. 



2007 



2009 



Rein et al.ll2010l ) . This is straightforward and 



precise, but t he cal culation is very time-consuming and the practical number of particles is limited. 
Tanga et al.l (|2004l ) performed A^-body simulations of a gravitationally unstable disk with a local- 
shearing box. They investigated the formation of clumps of planetesimals at t30 AU by the gravita- 
tional instability. They considered the drag force from the background gas. They showed that the 
planetesimal clumps form owing to the gravitational instability, which correspond to planetesimals 
in our calculation. In their simulation, the optical depth is smaller than that in our calculation and 
the particle disk is unstable even initially. We consider the formation of planet esimals at 1 AU 



thus we use the large optical depth, and assume the dust layer is initially stable. iMichikoshi et al 
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(|2007l ) performed a set of local simulations of the self-gravitating collisional particle disks without 
gas using iV-body simulation (hereaft er Paper I). Dust particle dynamics is calculated with the 
Hill equations in a local-shearing box ([Wisdom &: Tremaindll988l ). They adopted the rubble pile 
model (hard and soft sphere models) as collision models. They found that the for mation proc ess is 
divided into three stages: the formation of non-axisymmetric wake-like structures (|Salolll995l ). the 
creation of aggregates, and the rapid collisional growth of the aggregates. The m ass of the largest 
aggregates is larger than the mass predicted by the linear perturbation theory. iMichikoshi et al 



(j2009j) adopted the alternative model of collisions, 'accretion model' (Paper II). In the accretion 
model, the number of particles decreases as the calculation proceeds; thereby this model enables 
us to perform large-scale and long-term simulations. They obtain e d the final mass of planetesimals 
as a function of the size of the computational box. iRein et al.l (|2010l ) performed the numerical 
simulations to investigate the validity of the super-particle approximation. They considered the 
various physics, such as the gas drag, the self-gravity, the physical collision and the turbulence. 
They treated the collisions as the hard sphere model. They investigated the numerical requirements 
to study the gravitational instability. 

To understand the basic dynamics of gravitational instability, we neglected the effect of gas 
in Papers I and II. However, it is obvious that the gas plays an important role in planetesimal 
formation when dust particles are small. The interaction with gas through drag is dominant in 
the dynamical evolution of par ticles. Many autho r s investigated the effect of gas. The particle s 
drift radially due to gas drag (jAdachi et al.l Il976l : IWeidenschillingj Il977l ; iNakagawa et al.l Il986l ) . 
As the sedimentation of dust aggregates towards the midplan e proceeds, the vertical velocity 
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9j). The Kelvin-Helmholtz instability makes the dust laye r turbulent. The turbulent 



gas p revents the concentration of dust ([W eidenschill ing &; Cuzzil Il993l ) , and helps the concentra- 
tion (jBarge &; Sommerialll995l ; iFromang fc Nelsonl 120061) . The interaction between gas and dust 
causes the streaming instability (jYoudin fc Goodmanll2005l ). By the stre aming instability, the dust 



Johansen k, Youdin 2007; Johansen et al 
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grains concentrate strongly and the large Ceres- s ized p lanetesimals form (|Youdin h Johansenll2007l ; 

j). The loss o f angular momentum due to the 



200 



gas drag helps gravitational instability ( Ward 1976 : Youdin 2005al lbT). The time-scale for the dissi- 



pative instability is relatively slow, but the dust layer may be unstable even when the Toomre's Q 
value is larger than unity. As a first step toward understanding the effect of gas on gravitational 
instability, we introduce gas as a background laminar flow in the present paper. The aim of the 
present paper is to investigate the effect of the laminar gas on the planetesimal formation through 
gravitational instability. The laminar flow causes the dissipation of the kinetic energy and thus 
radial migration of dust. 

In §2, we summarize the results of the linear analysis of the gravitational instability under gas 
drag. In §3, we describe the simulation method and the initial condition. In §4, we present our 
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numerical results. §5 is devoted to discussions. We summarize the results in §6. 



2. Dispersion Relation for Gravitational Instability under Gas Drag 



The li near stability analyses of the dust layer with a finite thickness were performed by ISekiya 



(|l983l ) and lYamoto k, Sekiyal (120041 ) . They assumed that the size of dust aggregates is sufficiently 
small, in other words, dust is fully coupled with gas, thus they treated gas and dust as one fluid. 
That introduces no dissipation due to gas drag. They considered the vertical structure of the dust 



(1981 


), and 


Youdin 


(2005a 



section, we summarize the essence of their results. IWardl (|1976l ) and lYoudinl (|2005al ) used one 
component model, i.e., they neglected the dynamics of gas flow, but they consider the gas drag from 
the stationary background gas flow. Using their results, we discuss conditions for the gravitational 
instability under the influence of gas drag. 



2.1. Dispersion Relation 

We adopt the isothermal equation of state for the dust layer: 

P = c 2 p, (1) 

where p is the pressure, c is the isothermal sound speed, and p is the density of the dust layer. 
In this equation of state, the velocity dispersion of dust is always constant. Clearly, the velocity 
dispersion must not be constant. It changes owing to several processes such as gas drag, inelastic 
collisions, and gravitational scattering. This assumption gives us the simple analytical expression of 
the dispersion relation, with which we can grasp the nature of dissipative gravitational instability. 

We assume that the gas velocity is equal to the Kepler velocity. We consider the reference 
point that is at the semi-major axis do- The Kepler angular velocity of the reference point is 
= a/ GM s /cl'q where M s is the mass of the central star and G is the gravitational constant. We 
introduce the local Cartesian coordinates in which the x-axis is directed radially, the y-axis follows 
the direction of rotation, and the z-axis follows the direction perpendicular to mid-plane. We define 
v = {v x ,Vy,v z ) as the deviation velocity field from the local Kepler velocity. 

The basic equations for the dust layer in this frame are given by 

_^ + ( v • v)v = — £ + (2v y n, --v x n, - z n 2 ) + ~n x — - - — v - v<p, (3) 

at p 2 2 dy t stop 

V 2 = AirGp, (4) 
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where <j> is the gravitational potential due to the self-gravity. The simulation results in §4 show 
that the density fluctuation is not axisymmetric. However, for the sake of simplicity, restricting 
ourselves to the axisymmetric mode (d/dy = 0), we carry out a normal mode analysis in the form 
of exp(— i(ujt — kx)) where k is the wave number in the x-direction and uj is the growth rate. 



Integrating toward z, we obtain the equation of perturbed quantities Si, v x \, v y i, 4>\: 

— iuT,i + ikT,QV x \ = 0, 



iuiv x i = 2f2n 



ILOVyl 
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^stop 
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t 



stop 



-Vyl, 



27rGSi 



(5) 
(6) 

(7) 
(8) 



where £ is the surface density of the dust layer. lYoudinl (|2005al ) considered the finite thickness 
model using the softening term (|Vandervoortlll970l ) . Poisson equation with the softening parameter 
is 

- 2<irGY, 1 T(kh) = k<t>u (9) 
where h is the scale height of the dust layer and T(kh) is the softening parameter: 

1 



T(kh) 



1 + kh 



(10) 



For nontrivial solutions, the determinant of the coefficient matrix must vanish. Thus the 



following dispersion relation is obtained (jYoudinll2005al ): 



2 V + ( + -J- 



2itkG^oT(kh) + c 2 k 2 ) fi + 



^stop 



"-stop 



c 2 k 2 - 2TrkZ n G 



t 



stop 



(11) 



where [i = —iu. The mode for $l(p) > is unstable, where 3ft(X) is the real part of a complex 
number X. In the gas-free limit (t s top - > oo) and the thin disk limi t (kh — > 0), th e dispersion relation 



conve rges to the conventional dispersion relation of a thin disk (jToomre 



1964 



Goldreich &; Ward 



19731 ). Equation (jlip can be written in a non-dimensional form scaled by the length c/fl and the 
time Q^ 1 : 



- 3 + I 



1 



ti 



(12) 



<-stop y "stop 

where we used Toomre's Q given by Q = £Ic/(itGY<q). A tilde denotes non-dimensional quantity 
hereafter. 

A dissipative dust layer is secularly unstable for long- wavelength modes although Q > 1. The 
phys ical mechanism of the instability is explained as follows (boodman fc PindorhoOfi l IChiang &: Youdinl 
2010 ). We consider a thin axisymmetric overdense ring and the thick back-ground gas. From the 
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force balance of the self-gravity, we find that the dust rotates faster at the outer edge of the ring. 
Thus, the drag force at the outer edge is strong. This causes the inward drift. Conversely, the dust 
at the inner edge drifts outward. These effects shrink ring radially. 

Figure [U shows the dispersion relation of the dust layer under gas drag for Q = 0.5. We used 
the thin disk model (h = 0). The stopping time is 0.2, 1.0, 5.0, and oo (gas-free). All models have 
the maximum growth rate at k = 2. We will see that the growth rate has the maximum value at 
k = 1/Q in £ )2.21 As the stopping time lengthens, the growth rate shrinks. In the gas- free model, 
long-wavelength modes are stable, while in the gas model, they are unstable. 



2.2. Growth Rate of the Most Unstable Mode 

As discussed in Q2.1\ the dust layer is always unstable for long-wavelength modes. However, the 
timescale of the instability can be longer than the dynamical timescale, because this is the secular 
instability. Therefore, in order to understand what happens, we have to compare the timescale of 
the instability with those of other processes such as the sedimentation, the radial migration, and 
the growth of dust aggregates. 

Here, we assume the thin disk approximation (h = 0) in order to obtain the simple analytical 
formula. For k < 2/Q, fl is a real positive value. The growth rate for the most unstable mode \1 is 
a function of k, t s t op , and Q. From Equation (fl~2l) . the partial derivative of Jl with respect to k is 
given as 

Q± = 2(l/Q-k)(fi + l) 

dk 4/2/4to P + l/*Lp + ~k 2 + 1 - 2k /Q + 2,p ' 

Using Equation (|12|) . we can show that the denominator on the right hand side of Equation (|13p 
is positive when k < 2/Q. Hence, the jl has the maximum value at k = 1/Q for a fixed i s t p and 
Q. The substitution of k = 1/Q into Equation (|12p provides the equation of the maximum growth 

late /imax- 

~3 2 ~2 ( 1 1 \ ~ 1 

Amax + r Amax + I ~n hi— yr^ I Amax — ~ ~Fyi ~ \ ' 

tstop \ ''stop ^ / Istop^T 



To obtain the condition that the gravitational instability dominates over the other process 
(symbolically denoted by "x"), we equate Tx with //max where Tx is the characteristic timescale 
of a process X, ans solve the resultant equation with respect to Q to obtain Q value: 



Qcrit.X = Tx K 



^stop(^~X ^stop) (15) 



\ il tov f\ + (fx + ^stop ) 



When Q < Qcrit,x> gravitational instability is a dominant process compared to the process X. On 
the other hand, when Q > Q CT it,X) the timescale for the process X is shorter than the timescale of 
gravitational instability. Although the dust layer is unstable, the process X is dominant. 
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Fig. 1. — The dispersion relation of the dust layer under gas drag for Q = 0.5. We used the thin 
disk model (h = 0). The stopping time is 0.2 (dotted curve), 1.0 (short-dashed curve), and 2.0 
(dashed curve). The solid curve denotes the dispersion relation for the gas- free model. 
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We can derive the leading term of the power series expansion with respect to t R ton for E quation 
(13) and find the asymptotic solution (|Wardlll976l ; lYoudinll2005al : IChiang Youdinll2010l ). When 
the gas drag is strong t s t op <C 1, the growth rate for the most unstable mode is 



/'max 



■'stop 



(16) 



Conversely, when the gas drag force is weak t stop 3> 1 and Q > 1, the growth rate for the most 
unstable mode is 



(Q 2 



It 



stop 



When Q < 1, the growth rate for the most unstable mode converges to the finite value s/l/Q 2 — 1 
as t s top —> oo, which is the maximum growth rate for the gas-free gravitational instability. 



3. Method of Calculation 



3.1. Model 



The method of calculation is the same as those used in Papers I and II ex cept for including gas, 
which is based on the method of the local si mulation of planetary rings (e.g.. IWisdom &: Tremaine 
19881 : lRichardsonlll994l : baisaka & Idalll999h . 



We describe quantities in the non-dimensional form independent of the semi-major axis do, 
and the mass of the central star M s by scaling the time by J1q 1 , the length by the Hill radius 
han = fH, and the mass by /t 3 M s , where h = (2m p /3M s ) 1//3 and m p is the initial mass of particles 
(jPetit &: Henonlll986l : iNakazawa fc Idal ll988). The initial mass of particles m p is assumed to be 
ident ical. The non-dimensional Hill equations for the particle i are given by (e.g.. INakazawa fc Ida 
19881 ): 



d 2 Xi 

~w 

d 2 Vi 
dt 2 

d 2 Zj 
dt 2 



dt Y rf 3 

dXi ™i r~ 

Efhj ,, 



1 



t 



t 



stop,i 



stop,i 

dm 

dt 



dxi 
dt 



' g.r 



stop,i 



dzj 
dt 



(18) 
(19) 
(20) 



where 



77tj, and t s top,j are the position, mass, and the stopping time of the particle i, , n 



1 i 7 



is the distance between particles i and j, and (v g 



,) is the velocity field of gas. 



We neglect the dust back-reaction on gas for the sake of simplicity. We assume that the laminar 
gas flow, which has the sub-Keplerian velocity field: 



(Vg X , Vgy, Vg Z ) 



3 

0,- 2 x 



Vdi£,0 



(21) 
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where —3x/2 is the Keplerian shear velocity, v^if is the difference from the Keplerian velocity. The 
stopping time i s t p,i depends on the size of particle: i s top,i = t s t O p,0 ('fni/'fhp) ' where a is the power- 
law index, and i s top,o is the stopping time of the particle wit h m p . In this paper, we adopt a 
and 2/3. The index a = 2/3 corresponds to Stokes' law (e.g.. lLandau &: Lifshita 1 19591 ) . 







The calculation box has the periodic boundary condition (e.g., IWisdom fc Tremaind Il988l ), 
The size of the box is taken as a square AX m x ^4A m where A m is the most unstable wavelength of 
gravitational instability and A is the non-dimensional parameter. 

The system is characterized by two non-dimen sional parameters, th e optical depth and the 
ratio of the Hill radius to the diameter of a particle (jDaisaka &: Idal ll999): 



3S 



^"2r p ' 



(22) 



(23) 



where p p and r p are the internal density and initial radius of the particle, respectively. The realistic 
value for "solid" dust particles is ( ~ 105.28 for p p = 2g/cm 3 and ao = 1AU. However, we adopt 
£ ~ 2 in this paper owing to the computational limit (see Papers I and II). We set the optical depth 
r = 0.1, which is approximately equal to the realistic value r = 0.19 for = 10g/cm 3 
and p p = 2g/cm 3 . 



r p = 20cm, 



Since our choice of £ is smaller than the realistic value, a particle used in the simulation is not 
a realistic dust particle but a 'super-particle' that represents a group of many small particles that 
have the same position and velocity. Thus the restitution coefficient e corresponds to the rate of the 
dissipation due to collisions between super-particles. The low value of £ means that the physical 
size of a particle in our simulation is very large. This implies that the effect of the gravitational 
scattering is relatively weaker than that of collisions. We investigated the £ value dependence in 
paper I and II in the narrow range (£ = 1.5 — 3.0) and we confirmed that the formation process does 
not change in this range qualitatively. However, we show that the collisional growth is important in 
the final stage where \zeta controls the growth rate. Thus, we should perform the large £ value. In 
the future work, we will concentrate on the collisional growth and perform the iV-body simulations 
with a realistic £ value. 

In the Hill coordinates, the Keplerian orbit is determined by six parameters: the position of 
the guiding center, ecce ntricity e, inclination i, and two phases for epicyclic and vertical oscillations 
(jNakazawa Sz Idalll988l ). The initial eccentricity and inclinatio n of particles are ass umed to follow 
the Rayleigh distribution. We fix a/ (e 2 ) / (i 2 ) = 2 according to llda fc Makinol (|1992l ). 

The other parameters are uniformly distributed, avoiding overlapping. 

We adopt the hard-sphere and accretion models as a collision model (Pape rs I and II). In th e 
hard-sphere model, the penetration of particles is not possible (e.g., paper I, iRichardsonl 1 19941 ). 
When a collision occurs, the particle velocity changes instantly. The relative tangential velocity 
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is conserved and the magnitude of the relative normal velocity is redu ced by a factor of e . We 



use the constant e. We also use the accretion model (e.g., Paper II, iKokubo et al.l |2000| ). In 
this model, when two particles collide, if the binding condition is satisfied, the two particles merge 
into one particle . The binding condition is described by the Jacobi integral J < (e.g., Paper II, 



Nakazawa k Idalll988T ). 



In the hard-sphere model, since a planetesimal consists of many small particles, it can be 
distorted in shape and break up. On the other hand, in the accretion model, a planetesimal 
cannot be distorted or break up. In this respect, the hard-sphere model is more realistic than 
the accretion model. However the accretion model has an advantage. First, since the number 
of particles decreases as the calculation proceeds, this model enables us to reduce the calculation 
cost. Second, we can easily handle size-dependent drag coefficient. In general, the drag coefficient 
depends on the size of a particle. In the hard sphere model, it is complicated to handle the size- 
dependent drag coefficient. In the accretion model, we can adopt the size-dependent drag without 
any difficulty. We have shown that results in the accretion model are the same as those in the hard 
sphere model in the gas- free model (Paper II). 



3.2. Test Simulations 

We compare the accretion model with the hard sphere model in the laminar gas flow. We 
performed simulations with the same model parameters except for collision models (models 100H0, 
100AC0, and 100AS0) in order to check the validity of the accretion model. The initial stopping 
time of particles is t s top,o = 1-0. In the accretion models, we calculate the drag coefficient using 
the constant model (a = 0) and the Stokes' law model (a = 2/3). The left panel of Figure [2] shows 
the time evolution of Q. In the phase where Q decreases, there are no clear differences among all 
models. On the other hand, in the phase where Q increases, the discrepancies appear among them. 
Toomre's Q value in the accretion model with a = 2/3 is larger than those in the hard sphere 
model and the accretion model with a = 0. In the Stokes' law model, as particles grow, the drag 
coefficient becomes small and thus the dissipation becomes inefficient, which leads to the relatively 
large velocity dispersion. 

The right panel of Figure [2] shows the maximum mass of planetesimals for all models. In the 
early stage, particles or aggregates do not grow. There the difference in collision models is not 
important. The difference can be seen in the late stage. In the Stokes' law model a = 2/3, the 
drag coefficient changes as particles grow, which affects the evolution of the maximum mass. In 
the hard-sphere model, the drag coefficient is constant although aggregates grow, therefore the 
evolution in the hard-sphere model must be similar to that in the accretion model with a = 0. 
However, the difference appears at t/t^ = 10. The maximum mass of the accretion model with 
a = increases, but that of the hard-sphere model does not change. At i/tx = 10, there are only a 
few particles or aggregates. The maximum mass of planetesimals is sensitive to the initial condition 
and noise. But, at t/t^ = 20, the maximum mass converges to the same value. The final mass of 



the planetesimal depends on the size of the computational domain (jMichikoshi et al 
dispersion of the mass of the final state is relatively small. 



20091 ). The 



In the early stage where no planetesimals form, the time evolution is quite similar in all models. 
In the late stage, once planetesimals form, there are no remarkable differences among all models. 
If we focus on the the evolution of statistical value, such as the velocity dispersion, the number 
of planetesimals, we can investigate the gravitational instability and planetesimal formation using 
the accretion model. Although the time evolution of the maximum mass of the planetesimal in 
the late stage depends on the initial noise, the maximum mass finally converge to the same value. 
Therefore, we adopt the accretion model to investigate the gravitational instability from now on. 

4. Results 

We performed 14 simulations with different disk models. The model parameters are summa- 
rized in Tabled) We fixed the following parameters: the optical depth, r = 0.1, the restitution 
coefficient e = 0.01, the ratio of the Hill radius to the diameter £ = 2.0, the initial Toomre's Q value 
Qinit = 3, and the size of the computational domain A = 6. The basic dependence of planetesimal 
formation of these parameters was investigated in Papers I and II. 

4.1. Time evolution 

Figure [3] shows the typical evolution of the simulation (model 100AS0). At t/tK = 0, and 0.2, 
particles are distributed randomly and uniformly in the computational box. No gravitational insta- 
bility occurs, thus we can see no remarkable structures. At t/tK = 0.4, the large non-axisymmetric 
density structures appear. Then, the gravitational instability seems to start. In the model where 
(4top = 0.25), the structure is not clear. For the short stopping time models, the velocity dispersion 
is small and the critical wave length is short. Thus, we cannot observe the large density structure 
clearly. At t/tK = 0.6, particles begin to grow in the dense region. At t/tK = 0.8, and 1.0, particles 
continue to grow, and the number of small particles decreases rapidly. The small particles are 
absorbed by large particles. The rapid collisional growth of particles continues until the one large 
planetesimal finally forms. 

The formation process of planetesimals in the laminar gas model is similar to those in the 
gas- free model (papers I and II). The gas drag from the laminar gas does not qualitatively change 
formation process. This process is divided into three stages: the formation of density structures, 
the creation of planetesimal seeds, and their collisional growth. 



The line ar stability an alysis of gas-free gravitational instability shows that the disk is stable 
when Q > 1 (lToomrelll964! ). As discussed in §3, the dust layer in the laminar gas disk is secularly 
unstable although Q > 1. However, the growth time of the dissipative gravitational instability 
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Table 1. Simulation parameters 



Model 


t 

^stop,0 


Collision 


ot 




010AS0 


0.10 


Accretion 


2/3 





025AS0 


0.25 


Accretion 


2/3 





050AS0 


0.50 


Accretion 


2/3 





100AS0 


1.00 


Accretion 


2/3 





200AS0 


2.00 


Accretion 


2/3 





400AS0 


4.00 


Accretion 


2/3 





1000AS0 


10.00 


Accretion 


2/3 





INFA 


oo 


Accretion 







100AC0 


1.00 


Accretion 








100H0 


1.00 


Hard 








100AC1 


1.00 


Accretion 





10 


100AS1 


1.00 


Accretion 


2/3 


10 


100AC2 


1.00 


Accretion 





20 


100AS2 


1.00 


Accretion 


2/3 


20 



Note. — Parameters £ s top,o> cc, and are 
the stopping time of particles, the power-law 
index for the stopping time t s top,i oc rhf, and 
the velocity difference from Kepler velocity. 
We fixed other parameters: the optical depth 
r = 0.1, the restitution coefficient e = 0.01, 
the ratio of the Hill radius to the diameter 
C = 2.0, the initial Q value Qinit = 3, and the 
size of the computational domain A = 6. 
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is longer than the sedimentation in the initial stage. Thus, the disk shrinks vertically. As the 
velocity dispersion and the scale height decrease, the growth rate of the dissipative gravitational 
instability increases. Then, the gravitational instability becomes a dominant process, and the 
wake-like structures appear. These structures fragment into the seeds of planetesimals. The seeds 
grow rapidly owing to mutual collisions. Finally, almost all mass in the computational domain is 
absorbed by only one planetesimal in this calculation. The size of the final planetesimal depends 
on the size of the computational domain (paper II). 



4.2. Effect of Stopping Time 

We assume that the difference in the rotational velocity from the Kepler velocity is equal to 
zero (v dif = 0) (models INFA, 010AS0 - 1000AS0). In the model INFA, t stopfi = oo, i.e., the drag 
term is neglected. In other models, the stopping times are t s top,o = 0.1,0.25,0.5,1.0,2.0,4.0, and 
10.0. We adopt the Stokes' law model, a = 2/3. 



4-2.1. Toomre's Q value 

The top panel of Figure [4] shows the time evolution of Toomre's Q values. In all models, 
Toomre's Q value initially decreases. In the gas-free model (model INFA), the kinetic energy 
is dissipated only by inelastic collisions. The dissipation ti mescale due to inelastic col lisions is 



approximately estimated by t c ~ l/(rO(l — e 2 )) ~ lO.Ofi^ 1 (jGoldreich Tremaindll978l ). In the 
laminar gas models, the decrease of Toomre's Q value is caused by the dissipation due to both the 
gas drag and inelastic collisions. 

The top panel of Figure [5] shows the dependence of the decay time tdecay.Q on t s t O p,0! where 
the decay time of Toomre's Q value idecay,Q is defined by 

^decay,Q ^min, 

QTt E 7t ' ( 24 ) 

where Q m i n is the minimum Q value at t = t mm ,Q- The decay time tdecay.Q is proportional to 
£stop,o for i s top,o < 1- Toomre's Q value is proportional to the velocity dispersion a x . The inelastic 
collisions and the gas drag decrease the velocity dispersion. If the drag is sufficiently strong, the 
decay time scale of the velocity dispersion is on the order of the stopping time. The timescale of 
the decay of the velocity dispersion due to inelastic collisions is about t c ~ 10.0S1 - 1 . Thus, when 
£stop,o > 1 as t s top,o increases, we cannot neglect collisions. Therefore the decay timescale of the 
velocity dispersion is smaller than the stopping time for i s t O p,0 > 1- 

In all models, Toomre's Q value has the minimum value Qmin at imm,Q- In the gas-free model, 
the minimum Q value is about 2 (papers I, and II), and according to the linear theory, the critical 
value of a thin disk is Q = 1. This is because the decrease in the velocity dispersion due to inelastic 
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collisions is sufficiently slow. However, as shown in Figure 01 if the gas drag is strong, Q value 
becomes smaller than unity. As the stopping time shortens, the minimum Q value decreases. For 
example, the minimum Q value is 0.15 for i s top,o = 0.25. 

The axisymmetric structure starts forming when Q value is minimum. For example, in the 
model where i s top,o = 1-0, Q value is minimum at t/tj£ = 0.3, and we can see the non-axisymmetric 
structure at i/ijc = 0.35. Particles are scattered by the non-axisymmetric structures, and the 
velocity dispersion of particles increases. 

If the decrease of the Q value is more rapid than the gravitational instability, the dust layer 
cannot collapse although the dust layer is gravitationally unstable. As Q value decreases, the 
timescale of the gravitational instability becomes short, therefore, the critical Q value is determined 
by t s top,o = *Gii where toi is the timescale of the gravitational instability. Here, we assume t s top,i = 
£stop,o because particles do not grow much before gravitational instability. From Equation Q15]). we 
obtain the following condition: 



2t 2 

stop.O _ (25) 



^ + *stop,0 



The minimum Toomre's Q value as a function of t s top,o is shown in the bottom panel of Figure [5j 
Strictly speaking, Equation (|25p is valid only when the disk is thin. Nevertheless, the analytical 
expression agrees with the numerical result. 



[.2.2. Roche Density 



We define the scale height of a dust layer as the mean square root of z, h = \^2(zf) 1 ^ 2 . We 
can calculate the mean density of the dust layer from the scale height, p = T,/h. The density at 
which the self-gravity is equal to the tidal force is called the Roche density, which is defined by the 
following equation (e.g.. lYamoto &: Sekival 12004 ): 



9 M s 



(26) 



The dust layer may collapse if t he dust layer density exceeds the Roch e density. The linear analysis 
gives the more precise criterion ()Sekiyalll983l : lYamoto &; Sekival 120041 ). However, the criterion does 
not change very much. Thus, we introduce given by 



Qr 



PR 
P 



9 htt 2 



4tt GS ' 

which corresponds to the non-dimensional scale height. 



(27) 



The bottom panel of Figure S] shows the time evolution of Qr. The time evolution of Qr is 
similar to Toomre's Q. The value Qr has the minimum value /i m i n at t m in,h in all models. But, 
the time evolution of Qr is slower than that of Q value. This result indicates that the standard 
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hydrostatic relation h ~ \/2g z /£1 is not satisfied. Figure [6] shows the time evolution of the ratio 
h/o~ z , which is on the order of the sound crossing time across a scale height. In the gas-free model, 
the ratio h/a z is nearly constant in t/% < 1- The hydrostatic relation is satisfied. In a laminar 
gas model , the ratio h/a z is not constant. As the stopping time decreases, the variation of the 
ratio becomes steep. Especially, those for t < 1% is prominent. In the model t stop = 0.25, the 
ratio h/a z increases for t < 0.08% • The maximum value is about 3.7. The sound crossing time 
is h/a z / ~ 1.4 in the gas-free model. When t s top,o < 1> the dissipation is faster than the sound 
crossing. Therefore, the velocity dispersion decays before the scale height decays, and the ratio 
becomes larger than unity. In Figure we can see the oscillation. The period is about 0.4%, and 
the amplitude decreases gradually. The reason for the oscillation is not clear. This may be due to 
some kind of the relaxation process. 

We cannot use the Roche density as the collapse criterion for t stop > 1. For t s top,o > 4.0, Qr 
is always larger than unity. But the large particles form. For t s top.o = 2.0, Qr becomes unity 
at t/tK = 2, but the large particles start forming at t/% = 0-6- For t stop < 1, the time when 
Qr = 1 corresponds to the time when the large particles start forming. When we derive the 
Roche criterion, we simply compare the self-gravity with the tidal force. However, we showed that 
the non-axisymmetric structure develops. The azimuthal motion may be important to discuss the 
formation of the wake-like structure. Therefore, the Roche criterion is not applicable. 



The decay time of Qr value tdecay,h is defined by 



tdecay,h — tm'm,h , , j (28) 

"-init "-min 

where h m \ t is the initial scale height of a dust layer. Figure [7] shows tdecay.h as a function of t s top,o- 
The decay time of Qr has a minimum value at t s top,o = 0.5. For t s top,o < 0.5, the decay time 
shortens with the stopping time. For t s top,o < 0.5, the strong coupled particles with gas fall to 
the midplane at the terminal velocity. The settling time is proportional to t^opo- On ^ ne other 
hand, for t s top,o > 0.5, the decay time length ens with longer stopp ing time. In this regime, the 



particles oscillate around the equatorial plane (INakagawa et al.lll986l ). The decay time scale of the 
amplitude is proportional to the stopping time. 



4-2.3. Number of Planetesimal Seeds 

Figure [8] shows the time evolution of the number of planetesimal seeds per the area A^, N p \. 
We define a particle whose mass is larger than Mii near as a seed of planetesimals, where Mi; near is the 
planetesimal mass predicted by the linear theory 7rS(A m /2) 2 . In all models, N p \ has a maximum 
value at t/% = 0.5 — 1. The number N p \ increases rapidly in the early stage, and decreases in 
the late stage. In the early stage, gravitational instability occurs and the small density fluctuation 
grows and many planetesimal seeds form. In this stage, the number of seeds increases. They form 
through the gravitational instability. Once many planetesimal seeds form, the mutual collisions 
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among them start. The number of planetesimals decreases rapidly in the late stage. In the final 
state, almost all particles are absorbed by a few particles. The number of planetesimals in the final 
state depends on the size of the computational domain (paper II). 

The peak of N p \ depends on the stopping time t s top,o- Figure [9] shows the maximum number 
of the planetesimal iV max as a function of the stopping time i s top,o- The maximum number of plan- 
etesimals iV max has a peak at t s top,o — 0.5. As discussed in ^4.2.2t the scale height has the minimum 
value at t s top,o — 0.5. The characteristic length scale of gravitational instability corresponds to the 
scale height. As the characteristic length of gravitational instability shortens, the number of the 
planetesimals formed increases. 



4-2.4- Mass of Planetesimals 

Figure [10] shows the time evolution of the mass of the largest planetesimal. In all models, the 
largest planetesimal grows rapidly in t/tK < 5, and the growth is stalled at t/tx — 5 — 6. The 
growth in the gas-free model is slower than those in the gas models. As discussed in §4.2. H the 
gas drag enhances the dissipation of the kinetic energy of dust. As the drag becomes stronger, the 
gravitational instability occurs earlier and the growth becomes more rapid. 

The final mass of the largest planetesimal is 30 — 50Mii ncar . The clear dependence of the largest 
mass on the initial stopping time cannot be seen. As seeds grow, the stopping time lengthens. In 
the final state, we can neglect the effect of the drag. Therefore, the final mass and the number do 
not depend on the initial stopping time. Large planetesimals sweep small particles in the rotational 
direction. Therefore the final mass depends on the size of the computational domain (paper II). 



4.3. Effect of Rotational Velocity Difference and Drag Law 



The difference in rotational velocity between dust and gas causes the radial migration of dust. 
By neglecting th e self-gravity and the time derivatiyes o f Equation (|18p and (|19p . we obtain the 



steady solution: (lAdachi et al. 



1976 



Weidenschillind Il977l ) 



2t s top£dif 
+ 1' 



t 



(29) 



stop 



Figure [TT] shows the time evolution of the radial velocity v x of the largest particle. In these 
models, the initial stopping time is i s top,o = 1-0- The power-law index and the velocity difference are 
(a,v di f) = (0,10), (2/3,10), (0,20), and (2/3,20) (models 100AC1, 100AS1, 100AC2, and 100AS2). 
In all models, the radial velocity v x is negative at t/t^ = 1. This corresponds to the radially inward 
drift. The rotational velocity of gas is slower than the Kepler velocity; dust experiences a headwind, 
which causes its inward migration. Therefore their radial velocities become negative. Although a 
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is different, if i^if is same, the time evolutions of v x are similar in the early phase. This is because 
the planetesimals do not grow sufficiently in t < lix- If we estimate the terminal velocity from 
Equation ([29]) by the initial stopping time, they are the radial velocity v- m = 10 for Vda = 10 and 
u m = 20 for Udif = 20. In models for a = 0, the drag coefficient of a particle is independent of 
its size, thus the radial velocity converges to v x = 10 for v^a = 10 and v x = 20 for v^a = 20; the 
radial migration does not stop. On the other hand, In models for a = 2/3, as planetesimals grow, 
the drag coefficient becomes small; the radial migration stops finally. The radial velocity oscillates 
with the period of about tx- This oscillation is caused by the epicycle motion. The gas drag is so 
weak that the oscillation is not damped. 



5. Discussion 
5.1. Comparison with Linear Analyses 



Youdinl (|2005al ) performed linear analyses of the dissipative gravitational instability of a dust 
layer. They studied the effect of the stopping time, the thickness of the dust layer, and the 
velocity dispersion. We summarized the essence of their results in $2j The dust layer is always 
unstable for long wavelength modes. Thus, we should compare the timescale of the gravitational 
instability with those of other processes. Under the thin disk approximation (h = 0), we derived 
that simple condition that the gravitational instability dominates the other processes. Though the 
disk thickness cannot be neglected and the gravitational wakes is not axisymmetric, this condition 
corresponds to the occurrence of the gravitational wakes, as shown in To discuss the formation 
of the self-gravity wakes, Toomre's Q value is more important than the Roche criterion Qr. 

They applied the growth rate of the gravitational instability to the gravitational collapse in the 
turbulent disks. They assumed the simple turbulence model and studied the particle response to the 
turbulence flow. The protoplanetary disk can be turbulent due to magnetorotational instability or 
Kelvin-Helmholtz instability. They suggested that the standard hydrostatic relation h — V2a z /Q is 
not satisfied in the turbulent disk and the dust disk can be very thick. The thin disk approximation 
is clearly not valid. The criterion of the formation of the self-gravitational wakes may change. We 
plan to study the effect of the turbulence in the next paper. 



5.2. Comparison with Hydrodynamic Simulations 



Yamoto &; Sekival (|2006l ) performed the two-dimensional numerical simulation of the gravita- 



tional instability of a dust layer. They found that the dust layer becomes extremely thin if the 
stopping time is long, such as i s top > 0.1. This is because the dust settling is faster than the growth 
of the gravitational instability. However, our results show that the gravitational instability occurs 
although the stopping time is long. 
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There are so me differences between ou r simulations and theirs. We treated the dust component 
as particles. But, lYamoto Sekival (|2006l ) treated it as a pressure-free fluid. If the stopping time is 
short and the velocity dispersion is small, we can treat dust as a pressure-free fluid. The pressure- 
free fluid model becomes irrelevant for longer stopping times. For i st0 p = 0.1, we can use the 
pressure-free fluid, thus this should not cause the difference. In addition, the velocity dispersion 
sta bilizes the gravitational in stability generally. We neglected the dust back-reaction on the gas, 
but lYamoto Sekival (120061 ) adopted two-component model and solved the interaction between 
them. If the dust density is high and the stopping time is short, the dust back-reaction on the 
gas should be important. The resulting gas flow ma y change the criter i on. W e performed three- 
dimensional local calculation. On the other hand, lYamoto Sekival (|2006l ) assumed that the 
the dust layer to be axisymmetric with respect to the rotational axis. The gravitational instability 



forms the non-axisymm etric structures such as the gravitational wakes (e.g.. lMichikoshi et al 



2007 



Wakita fc Sekivall2008l ). The motion in the azimuthal direction is important to study the structure 
of the gravitationally unstable disk. 

To understand the dynamics of planetesimal formation in gas, we need to perform three- 
dimensional two-component simulations. We should consider the dust back-reaction on the gas 
and solve the hydrodynamic equations consistently. We will treat the dust component as particles 
because this method is applicable although the dust is loosely coupled to gas. We plan to study 
this point in the following papers of this series. 



6. Summary 

We performed local iV-body simulations of the planetesimal formation through gravitational 
instability. In papers I and II, we adopted the gas-free model. However the gas is not negligible 
and is important to the formation of planetesimals. As a first step to understand the effect of gas, 
we considered gas as the background laminar flow. We neglected the dust back-reaction on the gas 
for the sake of simplicity. Laminar flow causes the dissipation of the kinetic energy and thus radial 
migration of dust. 

In §2, we summarized the results of the linear analysis of the dissipative gravitational instability. 
To handle the dispersion relation analytically, we imposed the assumptions: the isothermal equation 
of state, the infinitely thin disk, and the axisymmetric perturbation. Simulations show that these 
assumptions are not valid, but they help us to understand the physical nature of the dissipative 
gravitational instability. We can use the growth rate to understand the simulation results. The 
long wavelength modes are always secularly unstable although Q > 1. Therefore, we must compare 
with the timescale of the gravitational instability with other processes. We provided the equation 
(|15p for the critical Q value. If Q is smaller than the critical value, the gravitational instability is 
a dominant process. 

The numerical simulations show that the formation process of planetesimals is the same as that 
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in the gas- free model. The formation process is divided into three stages qualitatively: the formation 
of wake-like structures, the creation of planetesimal seeds, and their collisional growth. By the linear 
stability analysis, the dust layer in the laminar gas disk is secularly unstable although Q > 1. The 
growth time of the dissipative gravitational instability is slower than the dust sedimentation and 
the decrease of the velocity dispersion when the initial Q value is sufficiently large. Thus, the disk 
shrinks vertically and Q decreases. As the velocity dispersion and the scale height decrease, the 
growth rate of the gravitational instability increases. Finally, the gravitational instability becomes 
dominant. Then wake-like structures are formed by the gravitational instability. These structures 
fragment into planetesimal seeds. Seeds grow rapidly owing to the mutual collisions. Finally, almost 
all mass in the computational domain is absorbed by only one planetesimal in the calculation. 

We investigated the dependence of results on the initial stopping time. The Q value decreases 
in the initial stage, and it reaches the minimum. The decay time scale of Q is proportional to 
the stopping time for i stoPj o < 1, which is the decay timescale of the velocity dispersion. The 
minimum Q is determined by the balance between the dissipation and gravitational instability. At 
the minimum Q, the wake-like structures are formed. The time evolution of the ratio of the dust 
layer density to the Roche density Qr is similar to that of Q. However the decay time scale of Qr 
is longer than that of Q. The time when Qr = 1 or Qr = Q_R im i n does not correspond to the onset 
of the gravitational instability, such as the formation of the wake-like structures. This indicates 
that the Roche criterion is not applicable in this system. We investigated the time evolution 
of the number of planetesimal seeds. We define a particle whose mass is larger than Mij near as a 
planetesimal seed, where Mii near is the planetesimal mass predicted by the linear theory 7rS(A m /2) 2 . 
The number of planetesimals has a maximum value. After the gravitational instability becomes 
a dominant process, planetesimal seeds form. In this stage, the number of planetesimal seeds 
increases. However, in the late stage, the number of seeds decreases by the mutual collisions among 
them. The maximum of the number of the planetesimal seeds depends on the initial stopping 
time. The maximum number of planetesimal seeds has a peak at i s top,o — 0.5. We are unable to 
identify a clear difference in the mass of the largest particle in the final state. The final mass of 



the planetesimal depends on the size of the computational domain (jMichikoshi et al.ll2009l ) 



We also investigated the models where the gas rotational velocity is slower than the Kepler 
velocity. We confirmed that the planetesimal is able to form in the sub-Keplerian gas disk. In 
the initial stage, particles experience a headwind, which causes its inward migration. If we adopt 
Stokes' gas drag law, as planetesimals grow, the drag coefficient becomes small; the radial migration 
finally slows. 

Almost all mass is absorbed by the largest planetesimal. In the gas-free models, we showed 
that the final mass of planetesimals depends on the size of computational domain. To investigate 
the final mass of planetesimals in detail, we should perform larger scale simulations. This remains 
to be discussed further. In addition, in this paper, we assumed gas to be laminar flow. However, the 
gas may be turbulent. If the turbulence is strong, the particles are stirred up and the gravitational 
instability may be prevented. The timescale for the dissipative gravitational instability is long 
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(Ward 


1976: 


Youdin 


2005a) 



and the timescale for the instability is long, the dust layer may disrupt owing to the dissipative 
gravitational instability over time in the turbulent disk. In the next paper, we will investigate the 
particle response to the turbulence and the gravitational instability in the turbulent disk. 

Numerical simulations were carried out on the MUV system at Center for Computational 
Astrophysics, National Astronomical Observatory Japan. This research was partially supported by 
MEXT (Ministry of Education, Culture, Sports, Science and Technology), Japan, the Grant-in- Aid 
for Scientific Research on Priority Areas, "Development of Extra-Solar Planetary Science," and the 
Special Coordination Fund for Promoting Science and Technology, "GRAPE-DR Project." S. I. is 
supported by Grants-in-Aid (16077202, and 18540238) from MEXT of Japan. 
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Fig. 2. — The time evolutions of Toomre's Q value (left panel) and the maximum mass M max of 
planetesimals (right panel) for hard-sphere model (solid curve), accretion model for a = (dashed 
curve), and accretion model for a = 2/3 (short dashed curve). 
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Fig. 3. — Spatial distribution of particles in the xy plane in in the model 100AS0 at t/tx = 
0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 4.0, and 6.0 from left to right and from top to bottom. Circles represent 
particles and their size is proportional to the physical size of particles. 
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Fig. 4. — The time evolutions of Toomre's Q value (top panel) and the normalized scale height 
(bottom panel) for gas- free model (solid curve), t s top,o = 4.0 (dashed curve), t s top,o = 2.0 (short 
dashed curve) , t s top,o = 1-0 (dotted curve) , t s top,o = 1-0 (dash-dotted curve), and t s top,o = 0.5 
(dot-short-dashed curve) (models INFA, 400AS0, 200AS0, 100AS0, 050AS0, and 025AS0). 
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Fig. 5. — The decay time tdecay,Q of the Toomre's Q value (top panel) and the minimum Q value 
Qmin (bottom panel) as a function of t s top,o- The cross denotes the result of Numerical simulations 
(models 010AS0, 025AS0, 050AS0, 100AS0, 200AS0, 400AS0, and 1000AS0). In the top panel, the 
solid line shows the line tdecay,Q = istop,o- I n the bottom panel, the solid line shows the analytical 
estimation of the minimum Toomre's Q value expressed by Equation (|25p . 
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Fig. 6. — The time evolutions of h/a for gas-free model (solid curve), t s top,o = 4.0 (dashed curve), 
4top,o = 2.0 (short dashed curve) , t s top,o = 1-0 (dotted curve) , t s top,o = 0.5 (dash-dotted curve), 
and t st0 p,o = 0.25 (dot- short- dashed curve) (models INFA, 400AS0, 200AS0, 100AS0, 050AS0, and 
025AS0). 
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Fig. 7. — The decay time of the scale height as a function of i s top,o- The cross denotes the 
result of Numerical simulations (models 010AS0, 025AS0, 050AS0, 100AS0, 200AS0, 400AS0, and 
1000AS0). 
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Fig. 8. — The time evolutions of the number of the planetesimals for gas-free model (solid curve), 
^stop,o = 4.0 (dashed curve), i s top,o = 2.0 (short dashed curve) , i s top,o = 1-0 (dotted curve) , 
£stop,o = 0.5 (dash-dotted curve), and t s top,o = 0.25 (dot- short- dashed curve) (models INFA, 400AS0, 
200AS0, 100AS0, 050AS0, and 025AS0). 
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Fig. 9. — The maximum number of the planetesimals as a function of i s top,o- The cross denotes 
the result of Numerical simulations (models 010AS0, 025AS0, 050AS0, 100AS0, 200AS0, 400AS0, 
and 1000AS0). 
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Fig. 10. — The time evolutions of the mass of the largest planetesimal for gas-free model (solid 
curve), t s top,o = 4-0 (dashed curve), t s top,o = 2.0 (short dashed curve) , i s top,o = 1-0 (dotted curve) , 
tstopfl = 0.5 (dash-dotted curve), and t s top,o = 0.25 (dot- short- dashed curve) (models INFA, 400AS0, 
200AS0, 100AS0, 050AS0, and 025AS0). 
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Fig. 11. — The time evolution of the radial velocity v x of the largest planetesimal for models 100AC1 
(solid curve), 100AS1 (dashed curve), 100AC2 (short dashed curve), and 100AS2 (dotted curve). 
Note that the planetesimal whose velocity is plotted is not identical, i.e., the largest planetesimal 
changes as the other planetesimal grows. Therefore the velocity changes abruptly then. 



